Przedmiotem analizy był zbiór zawierający pomiary trzyletnich śledzi oraz stanu środowiska, w którym żyją. Dane były zbierane na przestrzeni ostatnich 60 lat. Celem analizy było odkrycie przyczyny spadku długości śledzi na przestrzeni lat.
W ramach analizy dokonano podsumowania rozkładów zebranych atrybutów. Policzono także korelacje między zmiennymi. Okazało się, że najbardziej skorelowane są atrybuty reprezentujące liczbę różnych gatunków planktonu. Najmocniej skorelowana z długością śledzi jest natomiast temperatura przy powierzchni wody. W połączeniu z analizą nauczonego regresora prawdopodną wydaje się hipoteza, że atrybut ten w największym stopniu przyczynił się do obserwowanego spadku liczby śledzi.
library(corrplot)
## corrplot 0.84 loaded
library(mice)
## Loading required package: lattice
##
## Attaching package: 'mice'
## The following objects are masked from 'package:base':
##
## cbind, rbind
library(ggplot2)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(caret)
library(reshape2)
library(knitr)
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
W celu otrzymywania powtarzalnych wyników ustalono stałą wartość ziarna, która zostaje wykorzystywana przez algorytmy korzystające z pseudolosowości:
seed <- 420
set.seed(seed)
Wczytywanie danych odbywa się z użyciem funkcji read.csv. Przekazanie znaku zapytania jako argumentu na.strings powoduje, że wszystkie brakujące wartości zostaną poprawnie oznaczone jako NA.
data <- read.csv("sledzie.csv", na.strings = "?")
W celu uzupełnienia brakujących wartości wykorzystano bibliotekę mice. Skorzystano z metody cart, która stosuje drzewa regresyjne w celu predykcji brakujących wartości:
data <- complete(mice(data, meth='cart', seed = seed))
##
## iter imp variable
## 1 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 1 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 1 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 1 4 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 1 5 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 2 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 2 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 2 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 2 4 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 2 5 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 3 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 3 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 3 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 3 4 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 3 5 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 4 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 4 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 4 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 4 4 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 4 5 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 5 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 5 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 5 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 5 4 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 5 5 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
Rozmiar zbioru:
nrow(data)
## [1] 52582
Liczba atrybutów:
ncol(data)
## [1] 16
Przykładowe wartości (w tym przypadku pierwsze 6 pomiarów):
head(data)
## X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr
## 1 0 23.0 0.02778 0.27785 2.46875 21.67333 2.54787 26.35881 0.356 482831
## 2 1 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 3 2 25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 4 3 25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 5 4 24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 6 5 22.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## cumf totaln sst sal xmonth nao
## 1 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 2 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 3 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 4 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 5 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 6 0.3059879 267380.8 14.30693 35.51234 7 2.8
Z pomocą funkcji summary można podejrzeć minimum, maksimum, średnią oraz 1. i 3. kwartyl dla każdego atrybutu:
summary(data)
## X length cfin1 cfin2
## Min. : 0 Min. :19.0 Min. : 0.0000 Min. : 0.0000
## 1st Qu.:13145 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778
## Median :26291 Median :25.5 Median : 0.1111 Median : 0.7012
## Mean :26291 Mean :25.3 Mean : 0.4460 Mean : 2.0257
## 3rd Qu.:39436 3rd Qu.:26.5 3rd Qu.: 0.3333 3rd Qu.: 1.7936
## Max. :52581 Max. :32.5 Max. :37.6667 Max. :19.3958
## chel1 chel2 lcop1 lcop2
## Min. : 0.000 Min. : 5.238 Min. : 0.3074 Min. : 7.849
## 1st Qu.: 2.469 1st Qu.:13.427 1st Qu.: 2.5479 1st Qu.:17.808
## Median : 5.750 Median :21.673 Median : 7.0000 Median :24.859
## Mean :10.002 Mean :21.219 Mean : 12.8084 Mean :28.421
## 3rd Qu.:11.500 3rd Qu.:27.193 3rd Qu.: 21.2315 3rd Qu.:37.232
## Max. :75.000 Max. :57.706 Max. :115.5833 Max. :68.736
## fbar recr cumf totaln
## Min. :0.0680 Min. : 140515 Min. :0.06833 Min. : 144137
## 1st Qu.:0.2270 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068
## Median :0.3320 Median : 421391 Median :0.23191 Median : 539558
## Mean :0.3304 Mean : 520367 Mean :0.22981 Mean : 514973
## 3rd Qu.:0.4560 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351
## Max. :0.8490 Max. :1565890 Max. :0.39801 Max. :1015595
## sst sal xmonth nao
## Min. :12.77 Min. :35.40 Min. : 1.000 Min. :-4.89000
## 1st Qu.:13.60 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
## Median :13.86 Median :35.51 Median : 8.000 Median : 0.20000
## Mean :13.87 Mean :35.51 Mean : 7.258 Mean :-0.09236
## 3rd Qu.:14.16 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
## Max. :14.73 Max. :35.61 Max. :12.000 Max. : 5.08000
Wykresy pudełkowe dla każdego atrybutu (z pominięciem identyfikatora i miesiąca) prezentują się następująco:
data_m <- melt(select(data, -X, -xmonth))
## No id variables; using all as measure variables
ggplot(data = data_m, aes(x=variable, y=value)) + geom_boxplot() + facet_wrap( ~ variable, scales="free", ncol = 4)
W celu odkrycia najbardziej skorelowanych zmiennych przetransformowano macierz korelacji w macierz trójkątną oraz wypisano wszystkie znaczące korelacje (o współczynniku korelacji większym niż 0.5):
cor <- cor(data, use = "complete.obs")
corUpper <- cor
corUpper[lower.tri(cor, diag = TRUE)] <- 0
corHighest <- subset(as.data.frame(as.table(corUpper)), abs(Freq) > 0.5)
corHighest <- corHighest[order(-abs(corHighest$Freq)),]
colnames(corHighest) <- c("Atrybut 1", "Atrybut 2", "Współczynnik korelacji")
corHighest
## Atrybut 1 Atrybut 2 Wspólczynnik korelacji
## 101 chel1 lcop1 0.9562203
## 118 chel2 lcop2 0.8861733
## 169 fbar cumf 0.8165530
## 187 cumf totaln -0.7077921
## 116 cfin2 lcop2 0.6532125
## 247 lcop1 nao -0.5503011
## 253 sst nao 0.5122988
## 185 fbar totaln -0.5075128
## 245 chel1 nao -0.5058363
Wykres prezentujący graficznie korelacje między wszystkimi zmiennymi prezentuje się następująco:
corrplot(cor, type="upper")
Można zaobserwować wysoką korelację między dostępnością różnych rodzajów planktonu. Spowodowane jest to prawdopodobnie tym, że pomimo różnic międzygatunkowych każdy rodzaj planktonu potrzebuje mniej więcej zbliżonych warunków do życia i czynniki które wpływają na wzrost lub spadek liczby planktonu jednego rodzaju wpływają w analogiczny sposób na występowanie innych gatunków. Część korelacji jest stosunkowo oczywista i wynika z powiązania między atrybutami. Przykładowo zmienne fbar i cumf dotyczą tej samej wartości (natężenia połowów w regionie) i występuje między nimi bezpośrednia zależność. Atrybutem najmocniej skorelowanym z długością śledzi jest sst oznaczający temperaturę przy powierzchni wody.
Bardzo duża liczba pobranych wyników stawia wyzwanie przed czytelnym przedstawieniem zmiany rozmiaru śledzia w czasie. Pierwszy z pokazanych wykresów dokonuje uśrednienia 200 obserwacji występujących obok siebie:
data_len <- select(data, X, length)
n <- 200
data_aggr <- aggregate(data_len,list(rep(1:(nrow(data_len)%/%n+1),each=n,len=nrow(data_len))),mean)[-1]
ggplotly(ggplot(data_aggr, aes(x=X, y=length)) + geom_line())
Drugie podejście korzysta z obecnej w paczce ggplot metody geom_smooth, która wygładza wykres. W tym przypadku jako argument funkcji przekazywany jest cały zbiór:
ggplotly(ggplot(data, aes(x=X, y=length)) + geom_smooth())
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Oba wykresy pozwalają zauważyć, że rozmiar śledzi przez pierwsze około 17 tysięcy pomiarów, a następnie stopniowo spada.
Aby móc porównywać między sobą różne modele regresji zdefiniowano 5-krotną walidację krzyżową powtórzoną dwukrotnie. Dodatkowo przed uczeniem modelu wartości atrybutów są normalizowane poprzez odjęcie wartości średniej wewnątrz atrybutu oraz podzielenie wyniku przez odchylenie standardowe:
fitControl <- trainControl(method = "repeatedcv", number = 5, repeats = 2, preProcOptions = c('scale', 'center'))
Celem analizy jest zbadanie jaki wpływ na rozmiar śledzi ma środowisko, w którym żyją. Z tego powodu modele regresji trenowane są na danych pozbawionych atrybutów związanych z czasem, czyli “X” oraz “xmonth”:
dataLearning <- select(data, -X, -xmonth)
W ramach analizy skonstruowano modele regresji korzystając z wielu metod: * Regresja liniowa:
modelLm <- train(length ~ ., data = dataLearning, method='lm', trControl=fitControl)
modelLm
## Linear Regression
##
## 52582 samples
## 13 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 2 times)
## Summary of sample sizes: 42063, 42066, 42066, 42067, 42066, 42066, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 1.362265 0.3207204 1.08217
##
## Tuning parameter 'intercept' was held constant at a value of TRUE
modelLars <- train(length ~ ., data = dataLearning, method='lars', trControl=fitControl)
modelLars
## Least Angle Regression
##
## 52582 samples
## 13 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 2 times)
## Summary of sample sizes: 42064, 42066, 42065, 42066, 42067, 42066, ...
## Resampling results across tuning parameters:
##
## fraction RMSE Rsquared MAE
## 0.050 1.576921 0.2046328 1.273636
## 0.525 1.381235 0.3059592 1.100364
## 1.000 1.362400 0.3206643 1.082284
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was fraction = 1.
modelRidge <- train(length ~ ., data = dataLearning, method='ridge', trControl=fitControl)
modelRidge
## Ridge Regression
##
## 52582 samples
## 13 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 2 times)
## Summary of sample sizes: 42065, 42066, 42065, 42066, 42066, 42065, ...
## Resampling results across tuning parameters:
##
## lambda RMSE Rsquared MAE
## 0e+00 1.362488 0.3205524 1.082218
## 1e-04 1.362488 0.3205528 1.082224
## 1e-01 1.375899 0.3074362 1.095358
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 1e-04.
modelRf <- train(length ~ ., data = dataLearning, method='rf', ntree=25, importance=T, trControl=fitControl)
modelRf
## Random Forest
##
## 52582 samples
## 13 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 2 times)
## Summary of sample sizes: 42066, 42066, 42066, 42065, 42065, 42066, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 2 1.188379 0.4831498 0.9407136
## 7 1.194283 0.4780936 0.9447520
## 13 1.196616 0.4761352 0.9466780
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 2.
Najlepszą trafność uzyskano wykorzystując regresor Random Forest. Wyliczona ważność atrybutów prezentuje się następująco:
varImp(modelRf)
## rf variable importance
##
## Overall
## chel1 100.000
## lcop2 75.653
## sst 61.480
## cfin2 58.180
## recr 45.469
## totaln 41.738
## chel2 26.446
## nao 22.441
## fbar 22.142
## cumf 21.827
## sal 19.187
## lcop1 3.233
## cfin1 0.000
Przedstawione wyniki są znormalizowane, co oznacza, że najbardziej znaczący atrybut posiada wartość 100, a najmniej znaczący wartość 0. Atrybutem, który w największym stopniu odpowiada za zmianę długości śledzia jest atrybut cfin2 oznaczający dostępność planktonu Calanus finmarchicus gat. 2. Trzy kolejne atrybuty w kolejności znaczenia to totaln (łączna liczba ryb złowionych w ramach połowu), sst (temperatura przy powierzchni wody) oraz fbar (natężenie połowów w regionie). Wykresy dla tych atrybutów prezentują się następująco:
plots_cfin2 <- ggplot(data, aes(x=X, y=cfin2)) + geom_smooth()
plots_totaln <- ggplot(data, aes(x=X, y=totaln)) + geom_smooth()
plots_sst <- ggplot(data, aes(x=X, y=sst)) + geom_smooth()
plots_fbar <- ggplot(data, aes(x=X, y=fbar)) + geom_smooth()
grid.arrange(plots_cfin2, plots_totaln,plots_sst, plots_fbar, ncol=2)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Pierwszym wnioskiem, który można wysunąć jest podobieństwo kształtów wykresów dla atrybutów totaln oraz fbar. Atrybut fbar oznacza natężenie połowów totaln łączną liczbę śledzi. Wykresy potwierdzają intuicję, która podpowiada, że zwiększone połowy będą skutkować obniżeniem liczby śledzi zwłowionych w ramach pojedynczego połowu. Zjawisko to można zaobserwować np. podczas pierwszych 10 tysięcy obserwacji, gdy meljąca intensyfikacja połowów skutkowała dynamicznym wzrostem liczby złowionych jednorazowo ryb.
W przypadku ptzebiegu atrybutu cfn2 możemy zaobserwować przebieg zbliżony do przebiegu długości śledzi, lecz z lekkim przesunięciem. Atrybut ten oznacza występowanie planktonu konkretnego gatunku. Możliwym wyjaśnieniem obserwowanych zmian jest fakt, że rosnąca liczba ryb powoduje zmniejszenie ilości dostępnego planktonu, który stanowi dla nich pożywienie. Tłumaczy to początkowy spadek ilości planktonu, gdy liczba śledzi rośnie oraz późniejszy wzrost, gdy liczba ryb maleje. Tajemniczy pozostaje jednak gwałtowny spadek następujący mniej więcej po 35 tysiącach pomiarów. Atrybut sst charakteryzuje dokładnie odwrotny przebieg - początkowy spadek zmienia się w połowie pomiarów w znaczący wzrost. Analiza korelacji pokazała dodatkowo, że między atrybutem sst oraz długością śledzia występuje silna znacząca korelacja ujemna. Pozwala to sądzić, że przyczyną, która w dużym stopniu spowodowała spadek liczby śledzi może być wzrost temperatury przy powierzchni wody.